Draft version February 2, 2008 

Preprint typeset using I^T^jX style cmulatcapj v. 6/22/04 



ON FOREGROUND REMOVAL FROM THE WILKINSON MICROWAVE ANISOTROPY PROBE DATA BY 
AN INTERNAL LINEAR COMBINATION METHOD: LIMITATIONS AND IMPLICATIONS 

H. K. Eriksen 1 ' 2 - 3 

Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 
N-0315 Oslo, Norway 



A. J. Banday 

Max-Planck-Institut fiir Astrophysik, Karl-Schwarzschild-Str. 1, Postfach 1317, 
D-85741 Garching bei Miinchen, Germany 



K. M. GORSKI 3 

Jet Propulsion Laboratory, M/S 169/327, 4800 Oak Grove Drive, Pasadena CA 91109 
Warsaw University Observatory, Aleje Ujazdowskie 4, 00-478 Warszawa, Poland 

AND 

P. B. LlLJE 1 

Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 
N-0315 Oslo, Norway 
(Received March 1, 2004- Accepted May 20, 2004) 
Draft version February 2, 2008 

ABSTRACT 

We study the Internal Linear Combination (ILC) method presented by the Wilkinson Microwave 
Anisotropy Probe ( WMAP) science team, with the goal of determining whether it may be used for 
cosmological purposes, as a template-free alternative to existing foreground correction methods. We 
conclude that the method does have the potential to do just that, but great care must be taken both 
in implementation, and in a detailed understanding of limitations caused by residual foregrounds 
which can still affect cosmological results. As a first step we demonstrate how to compute the ILC 
weights both accurately and efficiently by means of Lagrange multipliers, and apply this method to 
the observed data to produce a new version of the ILC map. This map has 12% lower variance than 
the ILC map of the WMAP team, primarily due to less noise. Next we describe how to generate 
Monte Carlo simulations of the ILC map, and find that these agree well with the observed map on 
angular scales up to £ ~ 200, using a conservative sky cut. Finally we make two comments to the 
on-going debates concerning the large-scale properties of the WMAP data. First, we note that the 
Galactic south-eastern quadrant is associated with notably different ILC weights than the other three 
quadrants, possibly indicating a foreground related anisotropy. Second, we study the properties of the 
quadrupole and octopole (amplitude, alignment and planarity), and reproduce the previously reported 
results that the quadrupole and octopole are strongly aligned and that the octopole is moderately 
planar. Even more interestingly, we find that the £ — 5 mode is spherically symmetric at about 3tr, 
and that the I = 6 mode is planar at the 2cr level. However, we also assess the impact of residual 
foregrounds on these statistics, and find that the ILC map is not clean enough to allow for cosmological 
conclusions. Alternative methods must be developed to study these issues further. 
Subject headings: cosmic microwave background — cosmology: observations — methods: numerical 



1. INTRODUCTION 

The first-year release of the Wilkinson Microwave 
Anisotropy Probe (WMAP; Bennett et al. 2003a) data 
set has presented the cosmological community with an 
extraordinarily rich source of high-quality information, 
allowing the constraint of specific cosmological models 
and the parameters which define them to percentage ac- 
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curacy. 

Nevertheless, there remains an important goal beyond 
such a statistical assessment of the Cosmic Microwave 
Background (CMB) sky, namely to build an accurate 
image of the last-scattering surface which captures the 
detailed nature and morphology of our universe, and 
not simply some best- fit ensemble averaged view of it. 
Impediments to this program include instrumental noise 
and systematic artifacts, and foreground emission from 
local astrophysical objects. On a fundamental level non- 
cosmological foregrounds can easily compromise any con- 
clusion regarding primordial physics unless properly ac- 
counted for, while on a practical level they complicate 
both algorithms and analyses. Methods for either remov- 
ing, suppressing or at the very least constraining fore- 
grounds are therefore of great importance, and, indeed, 
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direct attacks on the raw data are very rarely justified. 
Practically any analysis must consider sky maps which 
have been processed in some way, either by explicit fore- 
ground corrections, or by introducing a sky cut. 

The importance of foreground removal has been rec- 
ognized in the community for a long time (e.g., Ban- 
day & Wolfendale 1991; Readhcad & Lawrence 1992; 
Brandt et al. 1994; Tegmark & Efstathiou 1996; Tegmark 
1998; Tegmark et al. 2000), as has the preferred method 
for discriminating against such contamination, namely 
multi-frequency observations. While the CMB itself con- 
tributes equally to all frequencies (as measured in ther- 
modynamic temperature units) due to the black body 
nature of its spectrum, foregrounds are typically strongly 
frequency-dependent. One may therefore distinguish 
between foregrounds and genuine CMB anisotropy by 
studying how the signal varies with frequency. How- 
ever, detailed subtraction of foregrounds has generally 
required one of two assumptions to be made: either that 
all of the physical components of the foreground emis- 
sion and their spectral behavior are known, or that accu- 
rate templates of all of the components arc available and 
that the appropriate spectral behavior can be determined 
by fitting the templates to the available multi-frequency 
data. Neither method is easily adapted to accommodate 
real spatial variations in the spectral behavior of the fore- 
grounds. 

The WMAP project appears to have systematic is- 
sues under control, whilst considerable effort has been 
expended on foreground issues, and uncertainties may 
still remain. Recent detections of non-Gaussianity (Chi- 
ang et al. 2003; Coles et al. 2004; Eriksen et al. 2004a; 
Naselsky et al. 2003; Park 2004; Vielva et al. 2004), the 
continuing debate over the low amplitude of large angular 
scale anisotropy (see e.g., Efstathiou 2004), and a possi- 
ble preferred direction or alignment contained therein (de 
Oliveira-Costa et al. 2004; Eriksen et al. 2004b) may yet 
be affected by improvements of our ability to remove non- 
cosmological foregrounds. The WMAP satellite observes 
the sky at five frequencies (23, 33, 41, 61 and 94 GHz), 
and using at least in part this information the WMAP 
team have applied three different methods for removing, 
constraining or describing the foregrounds (Bennett et 
al. 2003b). 

The first method is to produce a so-called internal lin- 
ear combination map (from now on denoted ILC), which 
assumes nothing about the particular frequency depen- 
dencies or morphologies of the foregrounds. Instead, a 
CMB map is reconstructed by co-adding the data at the 
five frequencies (now convolved to a common angular res- 
olution of 1 degree) with a set of weights that minimizes 
the final variance of the map. The details of the non- 
linear method adopted to derive these weights have not 
been described by the WMAP team. In order to accom- 
modate spectral variability of the foregrounds, the sky 
has been partitioned into 12 separate regions and the 
minimum variance criterion applied to each in turn to 
determine the weights. Discontinuities between regions 
have been minimized by smoothing of the weights at the 
boundaries. The resultant CMB map is visually satis- 
factory but has complex noise properties, and indeed the 
WMAP team explicitly warns against its use for cosmo- 
logical analysis. Nevertheless, the map has been sub- 
jected to such studies in the literature, and indeed the 



WMAP team do use the resultant map as an input to 
their second foreground removal technique. 

This involves the application of a Maximum Entropy 
Method (MEM) in order to construct a model of the 
foregrounds, component by component. The strength 
of this method in principle is its ability to reconstruct 
the synchrotron, free-free and dust emission and their 
detailed frequency dependence on a pixel-by-pixel ba- 
sis. However, the initial stage of the analysis must still 
utilize templates for these dominant foreground compo- 
nents, and also establishes priors on their spectral be- 
havior by using the WMAP data at the five frequencies 
after correction for a CMB component as determined by 
the WILC method above. As we will see later, ILC-likc 
methods in general still allow some leakage of foreground 
signal into the CMB reconstruction, and whether this re- 
sults in any feedback is difficult to determine. Again as 
a consequence of complex noise properties, the resultant 
map has not been considered useful for cosmological pur- 
poses. Instead, the WMAP team has used the results to 
interpret the nature of the foreground emission. In par- 
ticular, they identify a dust-correlated component in the 
lower frequency (23, 33 and 41 GHz) channels with a 
spectral index of (3 <~ —2.5 for a spectrum of the form 
v 13 . This is physically interpreted as a flat spectrum syn- 
chrotron component in regions of star formation near 
the Galactic plane, rather than to emission from spin- 
ning dust, which had become the preferred solution to 
this anomalous, low frequency dust correlated emission. 
The issue remains open, but recent reanalyses by La- 
gache (2003) and Finkbeiner (2004) find evidence for the 
latter interpretation. The origin of this controversy lies 
simply with the fact that the component separation as 
implemented by WMAP is allowed only to produce a 
combined synchrotron/spinning dust solution at each fre- 
quency, with no attempt made to separately disentangle 
these two components. 

The final method for foreground correction is perhaps 
the simplest of the three, and it is also the preferred 
method for generating cleaned maps suitable for cosmo- 
logical analysis. Rather than inherently exploiting the 
frequency information contained in the data, one sub- 
tracts external templates of the various physical com- 
ponents (i.e., maps produced by non-CMB observations 
made preferably at frequencies where a specific compo- 
nent dominates the emission) with coupling coefficients 
determined by cross-correlation with the observed maps. 
This avoids the noise amplification which occurs when 
one co-adds the WMAP data alone, and has the added 
benefit that the resultant maps have well-known noise 
properties, provided that the templates themselves do 
not contribute significantly to this. Difficulties associ- 
ated with the method include uncertainties in the de- 
tailed morphologies of the templates as scaled to the 
wavelengths of interest, and the propagation of errors in 
the coupling coefficients into the error budgets of specific 
analyses. 

It should be noted that Tegmark, de Oliveira-Costa, 
& Hamilton (2003) have applied a generalization of the 
ILC method to the WMAP data. The basic idea is to al- 
low the weighting of each map to be scale-dependent by 
performing the analysis in harmonic space, the assump- 
tion being that this allows any spatial variations in the 
spectral dependence of the foregrounds to be adequately 
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tracked. It is not clear to what extent real variations 
project onto the harmonic eigenmodes of the analysis. 
As with the ILC method, complex noise properties re- 
sult, and so it is unlikely that this method is suitable 
for high precision cosmological analyses. In what follows 
we will denote the map as TCM - the Tegmark et al. 
cleaned map. 

In this pa per a new look is taken at the ILC method 
presented bv lBennett et all l|2003b|) . with the main goal 
of determining whether a map derived in this manner 
can be suitable for cosmological purposes. Specifically, 
we derive a new ILC map based on Lagrange multipliers 
(in what follows to be referred to as the LILC - La- 
grange Internal Linear Combination - map), which has 
12% lower variance than the WMAP science team's ILC 
(hereafter referred to as WILC) map. We then generate 
Monte Carlo simulations of this map by adding white 
noise and foreground templates to CMB realizations, and 
process these through our ILC pipeline. This allows us to 
quantify the efficiency of the ILC method, and realistic 
foreground residual estimates may be established. 

In the final section we rep e at the large-scale analy- 
sis of Ide Oliveira-Costa et all (|2004) both for our new 
LILC map and for the simulations, to assess the impact 
of residual foregrounds on these statistics. However, we 
study not only the quadrupole and the octopole, but also 
consider the properties of the I — 4, 5 and 6 modes. In 
fact, we find that the properties of the latter two are at 
least as intriguing as those of the quadrupole and oc- 
topole: the I — 5 mode is highly spherically symmetric, 
and the £ — 6 mode is planar. 

2. METHOD 

The ILC method as defined bv lBennett et all (|2003bf) is 
based on a simple premise: suppose there are k observed 
CMB maps at different frequencies (but with identical 
beams), and the aim is to suppress foregrounds and noise 
as far as possible. Each of the k maps may be written 
(in thermodynamic temperature) on the form T{v^) = 
7cmb + T Iesidual (v k ), where T C mb and T leaidua x(y k ) are 
statistically independent. Therefore, if one now forms 
the linear combination 

fe 

T = Y d W i T{v i ), (1) 

i=l 

and requires that 

k 

X> = 1, (2) 

the resulting map may be written as 

k 

T = T C mb + WiT rosidua i(^). (3) 

Thus, the response to the CMB signal is always unity 
since it is independent of the frequency, and the k — 1 
free weights may be chosen to minimize the impact of 
the residuals. Assuming the CMB component is statisti- 
cally independent of the foregrounds and the noise, one 
convenient measure for this is simply the variance of T, 

Var(T) = Var(T CMB ) + Var ™*T rcsidua i(^)) • (4) 

n=i ' 



The internal linear combination method may now be de- 
fined succinctly in terms of Equations Q and [21 where 
the weights are determined by minimizing the variance 
in Equation 2] 

We compute the ILC weights by means of Lagrange 
multipliers. Our Lagrange multiplie r procedure is s imi- 
lar to the approach taken by iTeemark et al.l l)2003|) for 
computing the harmonic space weights from which their 
map is const ructed . A usefu l review of this method is 
also given bv lTegmarld l)1998j) . The variance of T is seen 
to be a quadratic form in the weights Wi, and its min- 
imization under the constraint given in Equation [21 is 
therefore most conveniently carried out by means of La- 
grange multipliers. As shown in Appendix A the linear 
system of equations to be solved can be written on the 
following form 
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where A is an arbitrary constant, w = (u>i, . . . , Wk) T are 
the ILC weights, and 

Cn ee (AT; AT}) = — - ]T(T*(p) - f*)(TJ(p) - T 3 ) 

* p— 1 

(6) 

is the map-to-map covariance matrix. The solutions to 
this system are the usual inverse covariance weights, 

= ¥^r- (7) 

Z^jk ^jk 

If the foreground properties vary strongly over the sky 
as a result of spatially dependent spectral indexes, then 
the ILC method may perform rather poorly. To remedy 
this, one may subdivide the sky into disjoint patches, 
and compute indepen dent set of weights for each patch. 
iBennett, et alJ ( 2003b) divided the full sky into twelve re- 
gions, eleven covering the non-uniform regions of Galac- 
tic plane, and the last one covering the Kp2 region plus 
the well-behaved parts of the Galactic plane. We will 
study this particular partitioning more closely in ^ 

Using such a partitioning, the minimization of the vari- 
ance in Equation 01 is carried out for each region sepa- 
rately, and the final step is therefore to construct one sin- 
gle full-sky map from those individual patches. In order 
to suppress boundary effects [Bennett et al.l l)2003b|) gen- 
erated a mask (i.e., a full-sky map consisting of 0's and 
l's) for each patch, and convolved these masks by a Gaus- 
sian beam of 90' FWHM. This final ILC map was then 
constructed by first generating one full-sky map from 
each ILC weight set, as described above, and then they 
co-added these maps pixel-by-pixel with weights given 
by the apodized masks. We adopt the same method for 
suppressing boundary effects without modifications. 

3. SIMULATIONS, CALIBRATION AND PERFORMANCE 

Most cosmological CMB analyses are based on Monte 
Carlo simulations, which in most cases is the only 
straightforward method of taking into account such real- 
world nuisances as non-uniformly distributed noise, non- 
Gaussian beam profiles and complex Galactic cuts. If 
the ILC cleaned map is to be used for such purposes, 
one must be able to construct a Monte Carlo ensemble 
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b) Average; full sky, noisy 




c) Average; 12 regions, noisy d) Standard deviation; 1 2 regions, noisy 

Fig. 1. — Plots showing the pixel-by-pixel average and standard deviation (lower right panel) of the difference maps taken between the 
reconstructed and the input CMB maps. Each plot is computed from 1000 simulations. Upper left: The full sky is treated as one single 
region, and no noise is added to the simulations. Upper right: Same as upper left, but WMAP specific noise is included. Lower left: The 
sky is partitioned into the 12 WMAP ILC regions, and noise is included. Lower right: The standard deviation of the difference maps for 
which the sky is split into 12 regions, and noise is included. 



that reproduces the detailed properties of the observed 
map. In this section, we first discuss how to produce 
such an ensemble, and then we take advantage of the 
simulations to study the properties of the ILC method 
itself. 

3.1. The simulation pipeline 

Monte Carlo simulation of the ILC map amounts sim- 
ply to producing a set of k base frequency maps with 
similar properties to the observed data, which are then 
processed through the ILC pipeline. The ILC pipeline 
may then in many respects be regarded simply as one 
among many statistics we apply to our maps - the crucial 
part is not the ILC pipeline in itself, but the construction 
of the base maps. The only difference from main-stream 
simulation is that we in this case add foregrounds to the 
simulations, rather than subtract them from the obser- 
vations. 

The simulation process may be written in the following 
algorithmic form 4 : 

1. Simulate one CMB component for each realization 
based on some power spectrum, and convolve this 
with the appropriate channel-specific beams. 

2. Add channel-specific foreground templates. 

3. Add a channel-specific noise realization. At this 
stage the simulation comprises 10 sky maps which 

4 Although summarized specifically for the WMAP processing, 
the method can clearly be generalized to other multi-frequency 
experiments. 



mimic the observed properties of the 10 WMAP 
channels at 5 frequency bands. 

4. For each channel-specific realization, deconvolve 
the beam, and convolve to a common resolution 
corresponding to a Gaussian beam of 1° FWHM. 

5. Compute an average map for each frequency. 

6. Apply the ILC pipeline. 

The only subtle point in this prescription is how to han- 
dle foregrounds. Ideally we would like to have a perfect 
full-sky, noiseless foreground template at each WMAP 
frequency and for each significant foreground (e.g., free- 
free, synchrotron and dust), but unfortunately, no such 
templates are available. We are therefore left with a 
choice between two options. 

First, we may use th e Fink b einer and H aslam 
templates dHaslam et alJ 119821 iFinkbeinerl 120041: 
Finkb einer et alJ 119991) for synchrotron, free-free and 
dust emission 5 , together wit h the channe l specifi c 
weights listed in Table 3 of iBennett et all l|2003b|) . 
The channel specific weights are estimated through 
direct fits to the observed data, and are therefore free 
of any assumptions about the spectral parameters. 
Moreover, this method includes the contribution from 
the anomalous dust-correlated foreground, without the 
necessity to resolve the spinning dust controversy. In 
practice, the weighted sum over the three templates 
approximates the correct amplitudes very well. 

5 Versions of these maps as used in the WMAP analysis are 
available at http://lambda.gsfc.nasa.gov 
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TABLE 1 

Efficiency as a function of region 



Region 


/* 


/// 


fd 


"CMB 
ff no j SG 


Full KpO region 


0.02 ± 0.06 


0.20 ±0.13 


0.54 ±0.06 


6.2 ±0.5 


Full Kp2 region 


-0.01 ± 0.04 


0.15 ±0.09 


0.54 ±0.06 


6.0 ±0.4 


Full Kp4 region 


-0.01 ± 0.03 


0.14 ±0.07 


0.55 ±0.07 


5.9 ±0.3 


WMAP ILC Kp2 


-0.04 ±0.01 


0.08 ±0.03 


0.55 ±0.06 


5.6 ±0.2 


Full sky 


-0.03 ± 0.01 


0.10 ±0.02 


0.46 ±0.03 


5.2 ±0.2 



Note. — The residual foreground levels and signal-to-noise ratios as a function of region size. f s - synchrotron fraction 
relative to the canonical contribution at 22.8 GHz; /// - the free-free fraction relative to 33.0 GHz; fd - the dust fraction 
relative to 93.5 GHz. The numbers are computed from sets of 1000 Monte Carlo simulations. 



On the other hand, at low Galactic latitudes the tem- 
plates approximate the real sky very poorly because 
of the complexity of the foreground emission and real 
spectral variations close to the Galactic plane, thus if 
a full-sky simulation is required, they should not be 
trusted. Nevertheless, for the purposes of calibration 
of our method, such inaccuracies are unlikely to affect 
our primary conclusions concerning the efficiency of the 
ILC method. Moreover, as we will demonstrate, in this 
implementation of the ILC method the inner Galactic 
region will always be strongly polluted by foregrounds, 
and should not be used for cosmological analysis. 

A second option is to use the MEM maps provided by 
the WMAP team. The advantage of this method is a 
much better approximation to the true sky emission at 
low Galactic latitudes. Unfortunately, these maps are 
intrinsically noisy, and one would therefore have to com- 
pensate for this when adding noise to the simulations. 
As a result, we adopt the simple template method in 
this paper, which results in simulations having accept- 
able power spectra in the high-latitude region. In fact, 
the simulations are visually acceptable even in the inner 
Galactic region, having features very much resembling 
those seen in the observed ILC map. 

3.2. Sensitivity and response to noise and region 
definitions 

While the ILC method itself is simple to define, it is 
less clear how accurately it allows the removal of Galactic 
foreground emission. To quantify this, we utilize our 
simulation set containing constant and known levels of 
these foregrounds, reconstruct the CMB sky estimate for 
each simulation via the ILC method, and compare this 
to the input CMB component. 

For the initial comparison, we consider the idealized 
case of a full-sky noise-less analysis, including only fore- 
grounds and CMB. The results from this exercise are 
shown in the upper left panel of Figure terms of the 
average residual map computed from 1000 simulations. 
In this case the ILC method does an excellent job of re- 
moving the foregrounds, as the residuals are less than 
10 /iK even in the central Galactic plane, about 0.01% of 
the K-band amplitude. The remaining residual is caused 
by the fact that it is possible to find a solution with 
slightly lower variance than even the true solution. 

The upper right panel shows the results from a simi- 
lar full-sky simulation, but for which Gaussian, channel- 
dependent noise has been added to each realization. The 
effect is striking, indeed, as both the Galactic plane and 



the North Galactic Spur are now clearly visible. The 
explanation lies of course in the definition of the ILC 
method - the ILC weights are defined to minimize the 
variance of the output map. In the noiseless case, this is 
an optimization only with respect to the foreground tem- 
plates; for three templates, no variations in the spectral 
indices, and four free weights to adjust, this can be per- 
formed to high precision. However, the problem becomes 
more complicated with the introduction of noise, since 
the minimum variance criterion then implies a trade-off 
between instrument or foreground noise. As is seen in 
Figure ^ a higher level of foregrounds in a relatively 
small region of the sky is preferred over increased noise 
over the full sky. 

Obviously even the clean, high-latitude regions of the 
sky become polluted by this higher level of foregrounds 
near the Galactic plane when treating the full sky as 
one region. In order to avoid such spreading one may 
therefore choose to divide the sky into separate patches, 
each with rather homogeneous foreground properties, as 
described above. While this procedure works very well 
in practice, there are certain problems that one should 
be aware of. 

In the lower two panels of Figure^we have plotted the 
average (lower left panel) and standard deviation (lower 
right panel) of the residual maps, when the sky is divided 
into the 12 regions defined by the WMAP team. Over- 
all the average map looks quite similar to the full-sky 
case, but there are a few important differences, namely 
that the inner galactic plane has a significantly smaller 
amplitude, and that the blue "halos" around it are less 
saturated. On the other hand, a few free-free regions are 
now visible, which were efficiently removed when treating 
the entire galactic plane as one region. 

However, the two most interesting points in this re- 
spect are to be found in the lower right plot, which shows 
the standard deviation of the difference maps. First, the 
scanning pattern of WMAP is clearly visible in the high- 
latitude region. This indicates that noise is more im- 
portant than foregrounds in this region, and therefore 
the ILC method prefers to minimize this, rather than 
for instance suppressing the North Galactic Spur, which 
is clearly visible in the average plot. Secondly, region 
number 12 (to the very right in the plot) is associated 
with a very large variance and so the estimated CMB 
signal is not only biased in this region, but for all prac- 
tical purposes unknown. In fact, in a number of noise- 
less simulations we have carried out the ILC weight ma- 
trix is singular in this region, indicating that there is 
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TABLE 2 
The WILC and LILC weights 



Region 


Map 


K-ba,nd 


K a,- band 


Q-band 


V-band 


W- band 


i 


WILC 


0.10876 


-0.68367 


-0.09579 


1.92141 


-0.25072 




LILC 


-0.19401 


0.14004 


07702 


0.61214 


36480 


2 


WILC 


0.10818 


-0.67987 


-0.09017 


1.96859 


-0.30674 




LILC 


-0 06280 


-0.14738 


-0 1 3982 


1.31073 


03927 


3 


WILC 


-0.04074 


-0.28682 


0.08476 


1.16221 


0.08061 




LILC 


-0.11470 


1 5098 


-0 38520 


1 1 6396 


18496 


■1 


WILC 


-0.01847 


-0.25533 


-0.02607 


0.83919 


.46068 




LILC 


-0 056^4 


-0.01464 


-0 31 223 


93407 


0.44934 


5 


WILC 


0.18610 


-0.77416 


-0.32352 


2.33978 


-0.42820 




LILC 


900QQ 




9789^ 




ZL^^O 


6 


WILC 


-0.02158 


-0.21880 


-0.08224 


0.84851 


0.47412 




LILC 


-0 10223 


21 569 


-0.51767 


.90277 


0.50144 


7 


WILC 


0.11790 


-0.67740 


-0.09117 


1.94830 


-0.29763 




LILC 


-U.UOUO ( 


n nnm ^ 

-U.UUU-LO 






n n^i f»Q 


8 


WILC 


0.12403 


-0.67639 


-0.09653 


1.74992 


-0.10103 






0. 104y4 


-0. 89602 


O.U/ /43 


2.013 ( / 


-0. 35952 


9 


WILC 


0.10500 


-0.68438 


-0.09847 


1.90588 


-0.22803 




LILC 


-0.04577 


-0.27660 


-0.02097 


1.28849 


0.05484 


10 


WILC 


0.16911 


-0.91455 


-0.01204 


2.64536 


-0.88788 




LILC 


0.19380 


-1.16103 


0.37899 


2.26627 


-0.67803 


11 


WILC 


0.21951 


-0.96567 


-0.18077 


2.38740 


-0.46046 




LILC 


0.22200 


-1.03357 


-0.09824 


2.34490 


-0.43509 


12 


WILC 


0.11101 


-0.67501 


0.05268 


1.59101 


-0.07970 




LILC 


-0.06397 


-0.00907 


-0.46855 


1.92500 


-0.38342 



Note. — Comparison of the official WMAP ILC weights and the Lagrange multiplier weights as derived in this paper. 



simply too little information present here to recover the 
CMB signal. When adding noise the matrix becomes 
non-singular, and the procedure does yield a result, but 
the reconstructed field is likely to be a very poor approx- 
imation to the underlying CMB field. The important les- 
son to be drawn from this is that the size of the patches 
must be large enough to provide adequate support for 
CMB reconstruction. Region number 12 is too small to 
do this, and should therefore either be merged into the 
large high-latitude region, or extended. 

3.3. Efficiency considerations 

By assuming a fixed spectral index for each of the 
important foregrounds it is possible to obtain reason- 
able estimates of the residual foreground level in the 
ILC map. Suppose each significant foreg round may be 
written in the form Tf{y) — {v/voYSo (jBennett et alJ 
2003Q), wh ere vq is an arbitrary reference frequency for 
the particular foreground, Sq is the true foreground dis- 
tribution on the sky at that frequency, (5 is the spectral 
index and T is measured in antenna temperature. Then 
the residual foreground contribution to the ILC map may 
be written on the form 



residual 



k 



— a(vi) 
vo 



So — f ■ S , (8) 



where a(y) is the conversion factor from antenna to ther- 
modynamic temperature. Thus, the parameter / is sim- 
ply the fraction of residual foregrounds of that particular 



type in the ILC map, relative to the chosen reference fre- 
quency. 

For the simulations, we know both the exact CMB 
component and the noise contributions, and so we can 
also compute the CMB signal to noise ratio by taking 
the ratio of the rms of the input CMB component to the 
noise rms. The latter quantity is computed as follows 

k 



E2 2 

z=l 



(9) 



where c^ oise i is the variance of the ith input noise real- 
ization, which we know. 

The efficiency of the ILC method may now be quan- 
tified by computing these parameters from the Monte 
Carlo simulations. Such results are summarized in Table 
^for five different high-latitude regions (including differ- 
ent amount of foregrounds). For each quantity we list the 
mean and standard deviation, as computed from a set of 
1000 simulations. Three foregrounds are included here, 
namely synchrotron (/3 S = —2.70, fo.s = 22.8 GHz), free- 
free (/3ff = —2.15, vo s = 33.0 GHz) and thermal dust 
(/3 d = +2.20, vo A = 93.5 GHz). 

The most interesting conclusions to be drawn from this 
table are the following: First, the ILC method performs 
quite well with respect to synchrotron emission, indepen- 
dently of the particular sky cut. Second, the more area 
is included in the analysis, the better it does for free-free 
emission, implying that the main support for informa- 
tion on free-free lies close to the Galactic plane, which is 
a reasonable result. 
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Third, the ILC method performs quite badly with re- 
spect to dust - the residual amount of dust in the simu- 
lations is approximately half that of the W-band, a point 
which must be well understood when using the ILC map 
for foreground studies. We will return to this issue in the 
next section. 

Finally, we see that the signal-to-noise ratio increases 
when excluding more of the Galactic plane. This is again 
a manifestation of the competition between noise and 
foreground reduction. When less foregrounds are in- 
cluded in the region of interest, relatively more emphasis 
is put on the noise. Thus, one can easily find, somewhat 
paradoxically, that by manually excluding foreground 
contaminated regions from the analysis, the final amount 
of residual foregrounds increases, simply because the area 
of choice does not carry enough information to calibrate 
the ILC weights properly, and therefore the ILC method 
preferentially eliminates noise rather than foregrounds. 

4. APPLICATION TO THE WMAP DATA 

Table El lists the ILC weights for each region and for 
each fre quency band, both as computed bv lBennett et al.1 
( 20033) and by the La grange multiplier method. Figure 
El shows our version of the ILC map. 

A comparison of the two weight sets in Table El shows 
clearly that the differences between the two methods are 
significant. The corresponding effect on the sky of these 
different weights is shown in Figure 01 a), where the dif- 
ference between our map and the WILC map is plot- 
ted. The most notable features include the large blue 
area around the Galactic bulge, presumably indicating 
the different abilities of the two methods to reject some 
large-scale foreground structures, and secondly the resid- 
ual small-scale structure most likely indicating the differ- 
ent noise properties of the two maps. 

In Figure |3{b) the difference between our map and the 
TCM is plotted (the TCM map was convolved to a com- 



mon resolution of 1° FWHM before computing the dif- 
ference). There are no noticeable small-scale structures 
uniformly distributed on the sky, indicating similar noise 
properties between the two maps. However, there are 
larger scale residual foreground features present. Some 
point-source-like residuals are also present, which are as- 
sociated with known WMAP sources. 

In order to assess the potential impact of point sources 
on our method, we computed weights for the Kp2 region, 
both including and excluding the 700 point sources re- 
solved by WMAP. The effect is negligible, at most a one 
or two percent modification of the weights. Nevertheless, 
this comparison does serve to remind us that there will 
likely be point source residuals in any ILC-derived CMB 
map. 

Another picture of the same comparisons is given in 
Figure 01 Here we have plotted the full-sky power spec- 
trum of the WILC map, the LILC map and the TCM, to- 
gether with the best-fit running index spectrum. Clearly, 
our map agrees very well with the TCM up to about 
£ = 200, but diverges at smaller scales, where the effect 
of the TCM's narrower W-band beam becomes clearer. 
The WILC map, however, departs from the other two 
already at £ w 30, a difference which is most naturally 
interpreted as resulting from different noise properties. 

We now compare the observed LILC power spectrum 
to simulated spectra. Figure [S] shows the power spec- 
trum of the observed ILC map together with 1 and 2a 
confidence band computed from 1000 simulations; the 
spectrum in the left panel is computed from the full sky, 
whereas the conservative KpO mask has been imposed 
in the right panel so that what is shown is actually a 
pseudo-spectrum. In the full-sky case, we see that the 
observed spectrum matches the simulations very well up 
to about £ « 100, but rises more rapidly from about 
£ > 150. When constrained to the KpO region, the ob- 
served spectrum follows the simulations all the way up to 





Multipole component, / 

Fig. 4. — Comparison of the full sky power spectra of the WILC 
map (red), the LILC map (blue) and the TCM (green). Notice the 
excellent agreement between the two latter spectra up to I xs 200, 
whereas the WILC spectrum departs from the other two already 
attfts 50. 



i — 200, after which a very small bias toward high values 
may be seen. Thus, the simulations seem to approxi- 
mate the real sky satisfactory on the KpO region, while 
they underestimate the level of residual foregrounds in 
the inner Galactic regions. 

The defining criterion of the ILC method is of course 
minimum variance, and the rms of the high-latitude re- 
gion of the LILC is 68 /xK, while the corresponding num- 
ber for the WILC is 72 fiK. In other words, our set of 
weights results in 12% lower variance, and is therefore 
better as far as the minimum variance definition is con- 
cerned. However, this does not necessarily mean that 
the level of residual foregrounds is smaller. In this, the 
contrary is true: by computing the residual fractions of 
each foreground in the high-latitude region as described 
in the previous section, we find that our map actually has 
slightly more foreground residuals than the WILC; the 
fractional residual foreground levels in the high-latitude 
WMAP ILC Kp2 region of the LILC map are [-0.069, 
-0.011, 0.736], while for the WILC map they are [-0.027,- 
0.017, 0.424]. 

As noted in the previous section, the amount of resid- 
ual dust is high in the ILC maps - the method is able 
to remove only half of the dust present in the W-band, 
where the dust is the dominant foreground. This re- 



sult is thu s in excellent ag r eemen t with the findings pre- 
sented by Naselsk y et alJ ((2003) , which concludes that 
the cleaned maps contain residual foregrounds which 
mainly originate from the W-band. 

4.1. Quadrant and hemisphere weights 

As pointed out earlier, one of the major weaknesses of 
the ILC method is its inability to handle spatial varia- 
tions in the spect ral indices of the foreg rounds. To rem- 
edy this weakness iBennett et alJ (2003b) divided the sky 
into 12 disjoint regions, and computed one set of weights 
for each region. Out of those 12 regions, 11 lie within 
the Kp2 Galactic plane, while the rest of sky was treated 
as one singl e region. In light of the asymmetries recently 
reported by Erikse n et al.l l|2004a|) . we have partitioned 
the high-latitude sky yet further, and subsequently com- 
puted weights for the Galactic hemispheres and quad- 
rants individually. The results from these computations 
are shown in Table |3J 

We first consider the quadrant numbers (quadrants 
are defined by the standard Galactic reference system.) 
While the NW, NE and SW quadrant numbers are 
approximately internally consistent, the SE quadrant 
stands out in the Q and V bands. Thus, these num- 
bers both support an d ask question of the findings of 
lEriksen et alJ |2004aj) . Certainly, the earlier results are 
supported in the sense that there is an asymmetry in the 
WMAP data, possibly marginally aligned from north- 
west to south-east. However, large differences in the 
weight coefficients would be interpreted most naturally 
in terms of variations of the noise and foreground proper- 
ties, in apparent contradictio n to the frequency inde pen- 
dence demonstrated both bv lEriksen et al.| l)2004a|) and 
lEriksen et alJ l)2004bfl . Further investigation is certainly 
warranted, but it may yet be that foregrounds could play 
a role in explaining the observed asymmetries. 

Unfortunately, it is difficult to assess the significance 
of the variations in Table |21 properly, but we can make 
a few rough estimates. We have generated 1000 sim- 
ulated realizations, and computed quadrant weights as 
described above for each of these. Then, for each re- 
alization we find the maximum absolute difference be- 
tween any two quadrants, for each frequency. The results 
from this exercise are summarized in Table in terms 
of the mean, the standard deviation and the maximum 
value found in the simulations. Note that these num- 
bers are only meant to give a rough idea of the shape 
of the distributions, and not for setting confidence levels 
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Multipole component, / Multipole component, / 

(a) (b) 

Fig. 5. — Comparison of the observed LILC power spectrum to the simulated spectrum. The gray bands correspond to 1 and 2cr 
confidence bands, estimated from 1000 Monte Carlo simulations. The left panel shows the spectrum computed over the full sky, and the 
right panel shows the pseudo-spectrum of the maps when the KpO mask is applied. 



TABLE 3 
Hemisphere and quadrant weights 



Region 


K band 


Ka band 


Q band 


V band 


W band 


Full Kp2 region 


-0.19401 


0.14004 


0.07702 


0.61214 


0.36480 


Northern hemisphere 


-0.20611 


0.14837 


0.13262 


0.55371 


0.37140 


Southern hemisphere 


-0.18015 


0.12169 


0.03213 


0.66930 


0.35703 


North-west quadrant 


-0.19451 


0.13659 


0.09579 


0.56725 


0.39489 


North-east quadrant 


-0.24447 


0.21397 


0.20529 


0.53331 


0.29190 


South-west quadrant 


-0.19393 


0.07268 


0.26102 


0.39229 


0.46793 


South-east quadrant 


-0.16324 


0.17738 


-0.23334 


1.01469 


0.20451 



Note. — Weights computed from Galactic hemispheres and quadrants outside the Kp2 mask. 



TABLE 4 

Maximum absolute quadrant weight differences 



Frequency 


Mean 


Std.dev. 


Max 


WMAP 


K-band 


0.064 


0.028 


0.181 


0.081 


Ka-band 


0.165 


0.073 


0.459 


0.141 


Q-band 


0.169 


0.074 


0.505 


0.494 


V-band 


0.157 


0.071 


0.502 


0.622 


W-band 


0.179 


0.082 


0.496 


0.263 



Note. — The distribution of maximum absolute weight 
differences between any two Galactic quadrants, as computed 
from 1000 simulations. The observed WMAP values are 



shown in the right-most column. 



- the distributions are non-Gaussian, and counting stan- 
dard deviations is therefore meaningless. Nevertheless, it 
is obvious that the quadrant differences observed in the 
true WMAP data are inconsistent with the adopted fore- 
ground model described by the combination of three tem- 
plates and fixed spectral indexes, and as proposed by the 
WMAP team. The weights of the south-east quadrant 
are radically different from those of the other three in the 



Q- and V-bands; the maximum difference found in the 
1000 simulations in the V-band is about 80% that of the 
observed data. Whether this indicates a real foreground- 
or noise-related problem in the south-eastern quadrant 
is not clear from this analysis, but it does question the 
validity of treating the entire high-latitude sky as one 
single region. 

The hemisphere results of Table |3 are by no means 
as decisive as the quadrant results, as the weights are 
more or less consistent with each other. However, this 
may very well be a coincidence; the internal variations 
between the north-west and north-east quadrants are 
much smaller than between the south-west and south- 
east quadrants, and yet, the two corresponding averages 
are rather similar. 

5. IMPLICATIONS FOR AND STABILITY OF THE 
LARGE-SCALE MODES 

In this Section we consider what implications our new 
LILC map have for the current debate concerning the 
peculiarities seen at the very largest scales, in particu- 
lar the questions of the seemingly low quadrupole, the 
planar octopole and the alignment between the two, and 
we establish the uncertainties connected to each of these 
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TABLE 5 

QUADRUPOLE AMPLITUDES 



Measurement 


5T|( M K 2 ) 


p- value 


Best-fit running index spectrum 


869.7 




Hinshaw et al. cut sky 


123.4 


0.018 


WMAP ILC map all sky 


195.1 


0.048 


Tegmark et al. 


201.6 


0.051 


Efstathiou - WILC map 


223.0 


0.063 


Efstathiou - TCM map 


250.0 


0.080 


Legendre ILC map 


350.6 


0.153 


Legendre ILC map (quadrants) 


345.1 


0.149 



Note. - Results from the measurements of the 
quadrupole amplitude. The third column lists the probability 
of finding a lower quadrupole than that of the corresponding 
map, given the theoretical model value shown in the first row. 



m easurements. For these s t udies we adopt the statistics 
of Ide Oliveira-Costa et all l)2004|) . and compute 1) the 
probability of finding a lower quadrupole moment than 
the observed one, given the best-fit WMAP spectrum, 
2) the probability of finding such a strong alignment be- 
tween the quadrupole and octopole and 3) the probability 
of finding such planar multipoles as seen in the WMAP 
maps. Here we briefly define the various statistics, 
and re fer the interested reader to lde Oliveira-Costa et al.1 
( 2004) for details on how each quantity actually is com- 
puted. 

5.1. Definitions of statistics 

The first statistic is simply the multipole amplitude 
ST 2 , which is defined in terms of a spherical harmonics 
expansion of the map, 



T(n) = V(ft) = J2 a imYim(n). 



I, in 



The multipole amplitude is then defined as 
1(1+1) 1 , „ 



ST? 



2tt 21 + 1 ^ 



\air. 



(10) 



(11) 



The next statistic is based on the possibility to define 
a preferred axis, ft/, for each multipole, namely that axis 
which maximizes the angular momentum dispersion, 



<V|(n-L) 2 |^) = 



i 

^] m 2 \a [m \ 2 . 

m= — l 



(12) 



The alignment between two modes is then measured by 
taking the dot product of the two preferred directions. 

The computation of this quantity is carried out by 
computing the spherical harmonic coefficients in some 
coordinate system, and then rotating these in harmonic 
space. Since the harmonic space rotation matrices are 
simple to compute, the complete maximization proce- 
dure becomes relatively inexpensive even for a high- 
resolution map with several million pixels. The details 
on computing these rotation matrices are described by 
Ide Oliveira-Costa et alJ lj2"004|) 6 . 

6 A complex conjugate of the harmonic coefficients may 
be necessary to obtain th e same results as reported by 



Ide Oliveira-Oosta et all I2004D . depending on which definition of 
the spherical harmonics one chooses. 
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Fig. 6. — The observed multipole amplitude plotted against 
the true, foreground-free amplitude. The observed WMAP LILC 
value is marked by a horizontal solid line, while the diagonal line is 
meant to guide the eye only; in the case of perfect reconstruction, 
all dots would lie along this line. Note that there are generally 
more dots above the dashed line than below it, indicating that the 
ILC reconstructed spectra are slightly biased toward high values. 
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Fig. 7. — The observed angular separation between the preferred 
quadrupole direction, 112, and the preferred octopole direction, n.3, 
plotted against the true, foreground- free separation. The color 
of each dot indicates the quadrupole amplitude of the given real- 
ization. For clarity we only plot those points which either lie in 
the 0—20% interval (blue squares) or in the 80—100% interval (red 
dots). The horizontal line indicates the quadrupole value for the 
LILC map. 
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TABLE 6 

Alignment of the quadrupole and the octopole 



Measurement 


h 


62 


h 


63 


Angle 


112 ■ 113 


WMAP ILC map all sky 


278° 


69° 


236° 


63° 


12° 


0.955 


Tegmark et al. 


258° 


59° 


238° 


62° 


10° 


0.984 


Legcndre ILC map 


247° 


62° 


233° 


63° 


7° 


0.993 


Legendre ILC map (quadrants) 


245° 


61° 


231° 


63° 


7° 


0.992 



Note. — Results from measurements of the position of the preferred directions of the quadrupole and octopole moments 
(denoted by Galactic longitude and latitude), and the alignment between these. The right-most column lists the probability of 
finding a weaker alignment between the quadrupole and octopole directions. 



The third quantity we consider is the degree of pla- 
narity of a given mode. Two different statistics for thi s 
purpose are defined bv Ide Oliveira-Costa et all (|2004j) . 
one which maximizes the angular momentum of the 
mode, and one which maximizes the fractional power 
that may be put into an azimuthal mode. We choose 
the latter, which may be written explicitly in the follow- 
ing form, 

|a ; _;| 2 + |a n | 2 . . 

t = max — j . (13) 

n V 1(7, I 2 

L-mi=-l \ u lm\ 

The maximization is performed over all pixels in the map. 
5.2. Results 

In Table the amplitudes of the quadrupole moments 
are tabulated for four different maps: the WILC map, 
the TCM, our LILC map, and finally the LILC map for 
which the Kp2 region is divided into quadrants. As we 
can see from the numbers in Table[S]the LILC quadrupole 
is significantly larger than those observed in the WILC 
map and the TCM map. In fact, according to our map, 
the CMB quadrupole is low only at the 1 to 7 level, 
or, in other words, it is in perfect agreement with the 
model. However, these measurements are associated with 
large uncertainties. It is true that there is no estimator 
induced va r iance in these measurements, as discussed by 
Efstathiou (2004), since we have access to the full sky, 
but we do know that the ILC method does not remove 
foregrounds perfectly in the presence of noise, and this 
obviously affects the large-scale modes. 

To assess the uncertainties in these measurements we 
once again take advantage of our simulations, for which 
we know both the CMB component and the recon- 
structed map, and compare the first few low-£ amplitudes 
for each realization. These results are shown in Figure 
E3 Each black dot in these plots indicates the true vs. 
the reconstructed amplitude for one Monte Carlo real- 
ization, and in the limit of perfect reconstruction, they 
should therefore all lie along the diagonal line. However, 
noise and residual foregrounds do have a significant ef- 
fect on these measurements, as seen by the considerable 
scatter in each panel. 

The solid horizontal lines indicate the LILC ampli- 
tudes, for which we obviously only know the recon- 
structed values. In the case of the much debated 
quadrupole amplitude, we see that the observed value 
of 351 /iK 2 may in fact originate from a cosmological 
quadrupole over the range ~ 130 /iK 2 to - 600 /iK 2 , 
and it is therefore difficult to assign a great deal of sig- 
nificance to this result. It is interesting to note that the 
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Fig. 8. — The observed t- value plotted against the true, 
foreground-free t-value. This parameter is defined as the fraction 
of the power attributable to the azimuthal component an to the 
total power C;, maximized over all reference frames. The symbols 
and colors have the same meanings as in Figure 171 



WILC map, which contains less residual dust than our 
map, also features a lower quadrupole. The most appro- 
priate conclusion to draw is that residual foregrounds can 
modify the quadrupole significantly, and it is important 
to propagate the uncertainties in foreground modeling 
into errors on such low order modes. 

Giv en this fact, the approach taken by lEfstathioil 
( 2004) might prove more reliable if the foreground un- 
certainties are dominated by residuals in the Galactic 
plane. In this work, the low-£ power spectra of the WILC 
and the TCM maps are estimated on a cut sky (based 
on the WMAP Kp2 mask). The most likely quadrupole 
amplitudes are found to be 223<uK 2 an d 250/iK 2 , respec- 
tively. An analysis bv lBielewicz et aTl l|2004|) utilizing a 
power equalization filter to reconstruct the low-£ multi- 
poles from the high-latitude signal yields similar results. 
Thus, a cut sky approach yields slightly higher values 
than the corresponding full-sky analysis does. 

We now turn to the question of alignment between the 
quadrupole and the octopole. The results from these 
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TABLE 7 
Planarity of a few multipoles 







£ = 3 




e = 5 


£ = 6 




Data set 


t 


P 


t 


p 


t 


P 


WMAP ILC map 


0.930 


0.124 


0.366 


0.999 


0.769 


0.031 


Tegmark et al. 


0.942 


0.096 


0.372 


0.998 


0.783 


0.024 


Legendre ILC map 


0.934 


0.114 


0.374 


0.998 


0.806 


0.015 


Legendre ILC map (quadrants) 


0.948 


0.081 


0.375 


0.998 


0.794 


0.019 



Note. — Results from measurements of the degree of planarity of the three multipoles, £ — 3, 5, 6. The left column in each 
section shows how much of the total power in the mode is attributable to the an component, as measured in a coordinate system 
in which the preferred direction is defined to be the z-axis. The right column shows the probability of finding a more planar 
multipole, as compared to an ensemble of 10 000 Gaussian simulations. 



measurements are summarized in Table El In this case 
we find that the alignment is actually stronger in the 
LILC map than in the WILC and TCM maps, having a 
probability as low as 0.7%. Again the associated variance 
is of great importance, and a scatter plot for these mea- 
surements are shown in Figure Each dot and square 
in this plot indicates the results from one simulation, for 
which we know both the input and output maps. In the 
limit of perfect reconstruction all dots should lie along 
the diagonal line. However, as seen from the very large 
scatter in this plot, it is clear that the ILC method does 
not reproduce the phases of the quadrupole and octopole 
modes accurately enough to justify a cosmological iden- 
tification. 

The colors in this plot indicate the quadrupole value 
of the reconstructed map for each realization, such that 
a realization marked by a blue square has an amplitude 
smaller than 80% of the simulations, and a realization 
marked by a red dot has an amplitude larger than 80% 
of the simulations. One would expect that an intrinsi- 
cally small quadrupole is more susceptible to foreground 
residuals than a strong one, and this is indeed the sit- 
uation. Further, as we have already seen, the observed 
LILC quadrupole value is very low indeed, and so the 
conclusion of the last paragraph is strengthened: An im- 
proved quadrupole estimate is required before we can at- 
tach cosmological significan ce to its propert i es. S imilar 
conclusions were rea ched in iBielewicz et al I l)2004|) and 
lHansen et a!~l (|2004lL 

Given that the LILC map contains more residual dust 
than the WILC and TCM maps and also features a 
stronger alignment between the quadrupole and the oc- 
topole, one may suspect that the alignment is driven by 
the dust component. However, no correlation was found 
between the alignment parameter t and the residual dust 
level /dust, or the two other foreground components. It 
is therefore difficult to conclude that the alignment is a 
direct result of residual foregrounds. 

Despite the arguments presented in the two previous 
sections, it is also worth noticing that there are very few 
dots below the LILC value even for the reconstructed 
maps, a fact which indicates that the ILC method does 
not seem to systematically introduce couplings between 
the quadrupole and octopole modes. Our results there- 
fore only demonstrate that there is a very large variance 
in this measurement, but not that there is a significant 
bias. The development of reliable cut-sky estimators of 
this feature seem to be of high importance. 



Finally, we turn to the issue of planarity and symmetry 
in the \ow-£ modes. In Table[7|we show results from mea- 
surements of the t-statist ic for the £ = 3, 5 and 6 modes . 
Interestingly, as noted by dc Olivcira-Costa et al. ( 2004), 
the octopole is planar roughly at the 1 to 10 level. The 
£ = 5 and 6 modes, however, are even more intriguing. 
The £ = 5 mode is highly spherically symmetric, and 
99.8% of the simulations have a larger t-value. On the 
other hand, the £ = 6 mode is strongly planar, with only 
1.4% of the simulations having a larger t- value. For com- 
pleteness, we note that the £ — 4 mode appears random 
in all respects in our analyses. 

In order to assess the foreground induced uncertainties 
in these measurements, we have plotted the observed t- 
value against the true value in figure |B1 It appears that 
the scatter dominates the results, and it is therefore dif- 
ficult to unambiguously conclude that the detections are 
truly cosmological in origin. However, we also see that 
the distributions are fairly symmetric about the diagonal 
line, indicating that the measurements are nearly unbi- 
ased. Residual foregrounds therefore seem to increase 
the variance in the measurements, but they do not ap- 
pear to introduce the sort of effects seen in the WMAP 
data. In this context, it is also important to notice that 
the observed WMAP values for £ = 5 and 6 are extreme 
compared to that observed in the processed ILC maps. 
Also, the power amplitude <5T| is very large, and conse- 
quently this mode should be quite robust against fore- 
ground perturbations. 

In figure|niwe have plotted the £ = 2, 3, 5 and 6 modes 
from the LILC map. Here we clearly see the origin of 
the effects discussed above: the planes determined by the 
peaks and troughs of the quadrupole and octopole appear 
to be very strongly aligned, while the degree of symmetry 
seen in £ — 5 is similarly striking. Finally, the £ = 6 mode 
is obviously highly planar, as seen by the very regular 
distribution of peaks and troughs. If such features can 
be unambiguously shown to be of cosmological origin, 
they may be indicative of new exotic physics. 

6. CONCLUSIONS 

The main goal of this paper was to study whether the 
ILC method is able to yield cosmologically useful maps, 
and if so, whether realistic simulations can be generated 
in reasonable time in order to calibrate the uncertainties 
associated with the use of such a map. The results pre- 
sented earlier suggest a cautiously positive conclusion - 
the ILC method does have the capability of producing 
relatively clean CMB maps without the use of external 





templates. Nevertheless, great care should be taken in 
the practical implementation of the method (e.g., the 
proper definition of the individual regions is a crucial 
step), and beyond this one needs to be highly aware of 
its limitations. 

On a more detailed level, we derived the equations 
for the ILC weights based on Lagrange multipliers, 
which were also discussed by Tegmark (1998). While 
a non-linear search algorithm is based on iterations, this 
method solves one single linear system of equations, and 
is therefore much faster. This is important when gen- 
erating Monte Carlo simulations. Subsequently, we dis- 
cussed how to produce realistic simulations of the ILC 
map, and used these simulations to study the proper- 
ties of the method itself, with particular emphasis on the 
sensitivity to noise and sky cuts. 

The method was applied to the real WMAP data, and 
the resultant LILC map was determined to have prop- 
erties similar to the TC M map, but somewhat different 
from the iBennett et all (|2003b) WILC map. We also 
computed ILC weights for four quadrants of the sky, and 
found that the south-eastern Galactic quadrant has sig- 
nificantly different properties than the other three, possi- 
bly she^dmgjiew light o n the asymmetry issue discussed 
bv lEriksen etall l|2004a|) . 

Finally, as a comment to the on-going debate on the na- 
ture of the large-angular scale anisotropy, we investigated 
the implications of the LILC map for estimates of the 
quadrupole and octopole modes, and found that the new 
quadrupole moment increases from 195 /iK to 351 /xK, 
which is a perfectly acceptable amplitude compared to 
the best-fit spectrum. However, the alignment between 
the quadrupole and the octopole is stronger in our map 



than in the WILC and the TCM. We also pointed out 
that the t = 5 and 6 modes are most peculiar in their 
symmetry properties, as only 0.2% of the simulations 
have a more spherically symmetric I = 5 mode than the 
WMAP data, and 1.4% a more planar £ — 6 mode. Fur- 
ther, since we have access to the full sky, these modes 
are all independent under the Gaussian, random-phase 
hypothesis, and the probabilities therefore accumulate 
quite straightforwardly. The major caveat, however, 
is that many of these measurements are derived from 
maps with complex foreground and noise properties, and 
definitive cosmological conclusions therefore remain elu- 
sive. Better foreground correction methods are required, 
or, alternatively, methods for studying the same prop- 
erties on a cut sky should be developed. This work is 
already under way, and will be published in a future pa- 
per. 

Returning to the ILC method, one may question 
whether the minimum variance criterion in itself is a 
meaningful measure of performance. As we have seen, 
this criterion implies a trade-off between suppressing 
noise and foregrounds, and moderate levels of fore- 
grounds are often accepted in order to suppress noise. 
For most practical cosmological analyses this is not likely 
to be acceptable; noise is more easily quantified than 
residual foregrounds. 

Note therefore that although we do provide a copy of 
the LILC map at H.K.E.'s home page 7 , we strongly ad- 
vise against using it for purposes beyond visual presen- 
tation, for which, of course, the official WILC map is 
perfectly acceptable. 

7 http:/ /www. astro. uio.no/~hke/cmbdata/WMAP_ILC_lagrange. 
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APPENDIX 

COMPUTING THE ILC WEIGHTS BY LAGRANGE MULTIPLIERS 

In this ap pendix we describe ho w to compute the ILC weights both efficiently and accurately by means of Lagrange 
multipliers. iBennett et al.l l)2fl03bfl do not specify how they carry out the minimization of Equation^in practice, other 
than stating that the minimum is found through a non-linear search. Thus, it is difficult to assess the accuracy of the 
final results they quote, as non-linear searches can often be plagued by convergence issues. However, again we remind 
the reader that the WMAP team only intended their ILC map to be used for visualization purposes, and obtaining 
high accuracy was therefore of little importance. Our goal, however, was to study whether this method may actually be 
used for cosmological purposes, as a credible alternative to the template correction method. In addition, Monte Carlo 
simulations was needed to fully account for the statistical noise properties of the method, therefore computational 
speed was a driving concern. 

Recall that the problem is to minimize the variance of a linear weighted sum over k frequency maps, as given in 
Equation under the constraint that the sum of weights equals unity. This latter constraint guarantees a correct 
response to the CMB component, while minimizing the variance suppresses residuals. 

Explicitly, the variance of the final map is given by 



Var(T) = (T 2 ) - (T) 2 
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where w = (wi, . . . ,Wk) T and 

Cy = (AT.AI)) = — E( T >) 

^ p—1 

is the map-to-map covariance matrix. 

Thus, the problem is simply to minimize a quadratic form, subject to the constraint given by Equation a task 
which is most conveniently solved by Lagrange multipliers. This problem can be restated slightly: First, we seek to 
minimize the following function, 

k 

/( w ) = E W iCii W 3 ( A3 ) 

under the constraint 

k 

3 (w) = 5> 2 = l. (A4) 

i=l 

In such cases the method of Lagrange multipliers tells us to look among those points, wo, which satisfies the following 
relation, 

V/(w ) = AV.g(w ), (A5) 

where A is an arbitrary constant. In other words, the extrema of /, subject to the constraint, g = 1, are just those 
points at which the gradients of / and g are parallel. 



T^T^ip) -T j ) 



(A2) 
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The partial derivatives of the function / are easily computed from Eauation lA3l and can be written on the following 
form, 

1 j=i 



The partial derivatives of g are obviously just unity. 
Thus, the extrema of / are found by simultaneously solving the system of k derivative equations given by Equation 
5] and the constraint in Equation [21 

(A7) 
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It can easily be shown that the weights solving this system are 

j—l ^ij 

u, i = — r — - 

and so we arrive at the usual inverse covariance weighting. 
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